******Parallel Trend Test*******
*** For Universities
use "D:\Data and Code\U_Direct", clear  

global xlist "size income lab age"
xtset id year, delta(1)
xtdescribe 

gen policy=year - 2019
tab policy
replace policy = -3 if policy <-3
replace policy = 3 if policy >3 

forvalues i = 3(-1)1{
	gen pre_`i' = (policy == -`i' & shidian == 1)
}

gen current = (policy == 0 & shidian == 1)

forvalue j = 1(1)3{
	gen post_`j' = (policy ==`j' & shidian == 1)
}

drop pre_1 //Drop the year before the policy


* Regression 1
xtreg h_sci pre_* current post_* $xlist i.year

* Figure 1
coefplot, baselevels ///
keep(pre_* current post_*) ///
vertical 
yline(0,lcolor(edkblue*0.8)) ///
xline(3, lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.75)) xlabel(,labsize(*0.75)) ///
ytitle("U: ESI Paper" , size(small)) ///
xtitle("", size(small)) ///
addplot(line @b @at) ///
ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
msymbol(circle_hollow) ///
scheme(s1mono)

* Regression 2
xtreg sci pre_* current post_* $xlist i.year

* Figure 2
coefplot, baselevels ///
keep(pre_* current post_*) ///
vertical 
yline(0,lcolor(edkblue*0.8)) ///
xline(3, lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.75)) xlabel(,labsize(*0.75)) ///
ytitle("U: Paper" , size(small)) ///
xtitle("", size(small)) ///
addplot(line @b @at) ///
ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
msymbol(circle_hollow) ///
scheme(s1mono)

*** For PRIs
use "D:\Data and Code\PRI_Direct", clear  

global xlist "size income lab age"
xtset id year, delta(1)
xtdescribe 

gen policy=year - 2019
tab policy
replace policy = -3 if policy <-3
replace policy = 3 if policy >3 

forvalues i = 3(-1)1{
	gen pre_`i' = (policy == -`i' & shidian == 1)
}

gen current = (policy == 0 & shidian == 1)

forvalue j = 1(1)3{
	gen post_`j' = (policy ==`j' & shidian == 1)
}

drop pre_1 //Drop the year before the policy


* Regression 3
xtreg h_sci pre_* current post_* $xlist i.year

* Figure 3
coefplot, baselevels ///
keep(pre_* current post_*) ///
vertical 
yline(0,lcolor(edkblue*0.8)) ///
xline(3, lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.75)) xlabel(,labsize(*0.75)) ///
ytitle("PRI: ESI Paper" , size(small)) ///
xtitle("", size(small)) ///
addplot(line @b @at) ///
ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
msymbol(circle_hollow) ///
scheme(s1mono)

* Regression 4
xtreg sci pre_* current post_* $xlist i.year

* Figure 4
coefplot, baselevels ///
keep(pre_* current post_*) ///
vertical 
yline(0,lcolor(edkblue*0.8)) ///
xline(3, lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.75)) xlabel(,labsize(*0.75)) ///
ytitle("PRI: Paper" , size(small)) ///
xtitle("", size(small)) ///
addplot(line @b @at) ///
ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
msymbol(circle_hollow) ///
scheme(s1mono)


******* Robust Test Using Alternernative Mathing Method ********
* PSM
use "D:\Data and Code\PRI&U_for_PSM", clear 
 global xlist "Size RD_Expenditure Lab Age center west"
 psmatch2 shidian $xlist, ///
 outcome(sci) radius logit caliper(0.05) ate ties common
 pstest $xlist, both graph
 
* DID
use "D:\Data and Code\PRI&U_for_Robust", clear  

global xlist "size RD_Expenditure lab age"
xtset id year, delta(1)
xtdescribe 

g time = (year >= 2019) & !missing(year)
g treated = (shidian >0)&!missing(shidian)
g did = time*treated

xtreg sci did  $xlist i.year, fe
xtreg h_sci did $xlist i.year, fe


******* Robust Test for Universities ********
use "D:\Data and Code\U_Direct", clear  

global xlist "size RD_Expenditure lab age"
xtset id year, delta(1)
xtdescribe 

g time = (year >= 2019) & !missing(year)
g treated = (shidian >0)&!missing(shidian)
g did = time*treated

*Alternernative Dependent Variable
g persci = sci/size
g perhsci = h_sci/size

xtreg persci did  $xlist i.year, fe
xtreg perhsci did $xlist i.year, fe

*Poisson Regression
tab year, gen(d_year)
tab id, gen(d_id)

xtpoisson sci did size RD_Expenditure lab d_year2-d_year5, fe
xtpoisson h_sci did size RD_Expenditure lab d_year2-d_year5, fe

*Winsorized data
winsor2 sci h_sci, replace cuts(1 99)

xtreg sci did  $xlist i.year, fe
xtreg h_sci did $xlist i.year, fe


******* Robust Test for PRIs ********
use "D:\Data and Code\PRI_Direct", clear  

global xlist "size RD_Expenditure lab age"
xtset id year, delta(1)
xtdescribe 

g time = (year >= 2019) & !missing(year)
g treated = (shidian >0)&!missing(shidian)
g did = time*treated

*Alternernative Dependent Variable
g persci = sci/size
g perhsci = h_sci/size

xtreg persci did  $xlist i.year, fe
xtreg perhsci did $xlist i.year, fe

*Poisson Regression
tab year, gen(d_year)
tab id, gen(d_id)

xtpoisson sci did size RD_Expenditure lab d_year2-d_year5, fe
xtpoisson h_sci did size RD_Expenditure lab d_year2-d_year5, fe

*Winsorized data
winsor2 sci h_sci, replace cuts(1 99)

xtreg sci did  $xlist i.year, fe
xtreg h_sci did $xlist i.year, fe

*************Placebo Test*************
***For PRIs
use "D:\Data and Code\PRI_Direct", clear

global xlist "size RD_Expenditure lab age"
xtset id year, delta(1)
xtdescribe 

*set matsize 10000
mat b = J(500,1,0) 
mat se = J(500,1,0) 
mat p = J(500,1,0) 

*repeated 500 times, regression on Paper
forvalues i=1/500{
	use "D:\Data and Code\PRI_Direct", clear
	xtset id year 
	keep if year==2016 
	sample 39, count   
	keep id  
	save match_id.dta, replace  
	merge 1:m id using "D:\Data and Code\PRI_Directa" 
	gen treat = (_merge == 3) 
	gen period = (year >= 2019) 
	gen did1 = treat*period
	reghdfe sci did1 $xlist ,absorb(id year) vce(cluster id)
		 
	* xtreg h_sci did1 $xlist i.year, fe
	 
	mat b[`i',1] = _b[did1] //* 系数矩阵
	mat se[`i',1] = _se[did1] //* 标准误矩阵
	mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did1]/_se[did1]))
}


svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .
label var pvalue1 pvalue
label var coef1 coefficient
keep coef1 se1 pvalue1  
save placebo_pri_sci.dta,replace   

*Plot
use placebo_pri_sci.dta,clear 
twoway(kdensity coef1,yaxis(1)) (scatter pvalue1 coef1, msymbol(smcircle_hollow) mcolor(blue),,yaxis(2)), ///
title("Placebo Test: PRI Paper" ) ///
xline(9.889,lwidth(vthin) lp(shortdash)) xtitle("Coefficients") ///
yline(0.0095, lwidth(vthin) lp(dash)) ytitle(kdensity of estimates) ///
legend(label(1 "kdensity of estimates") label(2 "p value")) ///
plotregion(style(none)) graphregion(color(white))


*set matsize 10000
mat b = J(500,1,0) 
mat se = J(500,1,0) 
mat p = J(500,1,0) 

*repeated 500 times, regression on ESI Paper
forvalues i=1/500{
	use "D:\Data and Code\PRI_Direct", clear
	xtset id year 
	keep if year==2016 
	sample 39, count   
	keep id  
	save match_id.dta, replace  
	merge 1:m id using "D:\Data and Code\PRI_Directa" 
	gen treat = (_merge == 3) 
	gen period = (year >= 2019) 
	gen did1 = treat*period
	reghdfe h_sci did1 $xlist ,absorb(id year) vce(cluster id)
	 
	* xtreg h_sci did1 $xlist i.year, fe
	 
	mat b[`i',1] = _b[did1] //* 系数矩阵
	mat se[`i',1] = _se[did1] //* 标准误矩阵
	mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did1]/_se[did1]))
}


svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .
label var pvalue1 pvalue
label var coef1 coefficient
keep coef1 se1 pvalue1  
save placebo_pri_hsci.dta,replace   

*Plot
use placebo_pri_hsci.dta,clear 
twoway(kdensity coef1,yaxis(1)) (scatter pvalue1 coef1, msymbol(smcircle_hollow) mcolor(blue),,yaxis(2)), ///
title("Placebo Test: PRI ESI Paper" ) ///
xline(9.889,lwidth(vthin) lp(shortdash)) xtitle("Coefficients") ///
yline(0.0095, lwidth(vthin) lp(dash)) ytitle(kdensity of estimates) ///
legend(label(1 "kdensity of estimates") label(2 "p value")) ///
plotregion(style(none)) graphregion(color(white))


***For Universities
use "D:\Data and Code\U_Direct", clear

global xlist "size RD_Expenditure lab age"
xtset id year, delta(1)
xtdescribe 

*set matsize 10000
mat b = J(500,1,0) 
mat se = J(500,1,0) 
mat p = J(500,1,0) 

*repeated 500 times, regression on Paper
forvalues i=1/500{
	use "D:\Data and Code\U_Direct", clear
	xtset id year 
	keep if year==2016 
	sample 39, count   
	keep id  
	save match_id.dta, replace  
	merge 1:m id using "D:\Data and Code\U_Directa" 
	gen treat = (_merge == 3) 
	gen period = (year >= 2019) 
	gen did1 = treat*period
	reghdfe sci did1 $xlist ,absorb(id year) vce(cluster id)
		 
	* xtreg h_sci did1 $xlist i.year, fe
	 
	mat b[`i',1] = _b[did1] //* 系数矩阵
	mat se[`i',1] = _se[did1] //* 标准误矩阵
	mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did1]/_se[did1]))
}


svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .
label var pvalue1 pvalue
label var coef1 coefficient
keep coef1 se1 pvalue1  
save placebo_u_sci.dta,replace   

*Plot
use placebo_u_sci.dta,clear 
twoway(kdensity coef1,yaxis(1)) (scatter pvalue1 coef1, msymbol(smcircle_hollow) mcolor(blue),,yaxis(2)), ///
title("Placebo Test: U Paper" ) ///
xline(9.889,lwidth(vthin) lp(shortdash)) xtitle("Coefficients") ///
yline(0.0095, lwidth(vthin) lp(dash)) ytitle(kdensity of estimates) ///
legend(label(1 "kdensity of estimates") label(2 "p value")) ///
plotregion(style(none)) graphregion(color(white))


*set matsize 10000
mat b = J(500,1,0) 
mat se = J(500,1,0) 
mat p = J(500,1,0) 

*repeated 500 times, regression on ESI Paper
forvalues i=1/500{
	use "D:\Data and Code\U_Direct", clear
	xtset id year 
	keep if year==2016 
	sample 39, count   
	keep id  
	save match_id.dta, replace  
	merge 1:m id using "D:\Data and Code\U_Directa" 
	gen treat = (_merge == 3) 
	gen period = (year >= 2019) 
	gen did1 = treat*period
	reghdfe h_sci did1 $xlist ,absorb(id year) vce(cluster id)
	 
	* xtreg h_sci did1 $xlist i.year, fe
	 
	mat b[`i',1] = _b[did1] //* 系数矩阵
	mat se[`i',1] = _se[did1] //* 标准误矩阵
	mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did1]/_se[did1]))
}


svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .
label var pvalue1 pvalue
label var coef1 coefficient
keep coef1 se1 pvalue1  
save placebo_u_hsci.dta,replace   

*Plot
use placebo_u_hsci.dta,clear 
twoway(kdensity coef1,yaxis(1)) (scatter pvalue1 coef1, msymbol(smcircle_hollow) mcolor(blue),,yaxis(2)), ///
title("Placebo Test: U ESI Paper" ) ///
xline(9.889,lwidth(vthin) lp(shortdash)) xtitle("Coefficients") ///
yline(0.0095, lwidth(vthin) lp(dash)) ytitle(kdensity of estimates) ///
legend(label(1 "kdensity of estimates") label(2 "p value")) ///
plotregion(style(none)) graphregion(color(white))






























                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              



